#ifndef NURBS_H
#define NURBS_H

#include "Curves/Interfaces/icurve.h"

template< class IPoint >
class NURBS : public ICurve< IPoint >
{

protected:
    /* Fields. */
    int degree;
    int num_of_knots;
    vector< double > knots;
    vector< double > weights;

public:
    /* Constructors. */
    NURBS( int degree );
    /* Destructors. */
    ~NURBS( );

    /* Methods. */
    virtual void AddPoint( const IPoint& point, double weight );
    virtual void AddPoint( const IPoint& point, int point_index, double weight );

    virtual void RemovePoint( );
    virtual void RemovePoint( int point_index );

    virtual IPoint Eval( double t ) const;
    virtual vector< IPoint > ComputePoints( int n ) const;

protected:
    void UpdateKnots( );
    int FindIndex( double t ) const;

};

/* Constructor. */
template< class IPoint >
NURBS< IPoint >::NURBS( int degree )
    : degree( degree )
{
    num_of_knots = 1;

    for ( int i = 0; i <= degree; i++ )
    {
        knots.push_back( 0.0 );
        knots.push_back( 1.0 );
    }

    UpdateKnots( );

}

/* Destructor. */
template< class IPoint >
NURBS< IPoint >::~NURBS( )
{
}

/* Methods. */
template< class IPoint >
void NURBS< IPoint >::AddPoint( const IPoint &point, double weight )
{
    ICurve< IPoint >::AddPoint( point );
    weights.push_back( weight );

    if ( this -> control_points.size( ) > static_cast< unsigned int >( degree ) )
    {
        /* Powiekszamy rozmiar listy wezlow. */
        knots.reserve( knots.size( ) + 1 );
        //knots.push_back( 1 );
        num_of_knots++;

        /* Aktualizujemy wartosci na skrajach przedzialow. */
        UpdateKnots( );
    }
}

template< class IPoint >
void NURBS< IPoint >::AddPoint( const IPoint &point, int point_index, double weight )
{

}

template< class IPoint >
void NURBS< IPoint >::RemovePoint( )
{

}

template< class IPoint >
void NURBS< IPoint >::RemovePoint( int point_index )
{

}

template< class IPoint >
IPoint NURBS< IPoint >::Eval( double t ) const
{
    assert( 0 <= t && t <= 1 );

    int j = FindIndex( t );

    /* t nie moze wpadac do tego przedzialu && wartosc nie wpada do zadnego przedzialu. */
    assert( static_cast< unsigned int >( j ) <= this -> control_points.size( ) - degree && j != -1 );

    IPoint **W;
    double **gamma, **alpha;

    W     = new  IPoint*[j + degree + 1];
    gamma = new  double*[j + degree + 1];
    alpha = new  double*[j + degree + 1];

    for ( int i = 0; i < j + degree + 1; i++ )
    {
        W[i]     = new  IPoint[j + degree + 1];
        gamma[i] = new  double[j + degree + 1];
        alpha[i] = new  double[j + degree + 1];
    }

    for ( int i = j; i <= j + degree; i++ )
    {
        W[i][0]     = this -> control_points[i];
        gamma[i][0] = weights[i];
    }

    for ( int k = 1; k <= degree; k++ )
    {
        for ( int i = j + k; i <= j + degree; i++ )
        {
            alpha[i][k] = (t - knots[i]) / (knots[i + degree + 1 - k] - knots[i]);
            gamma[i][k] = alpha[i][k] * gamma[i][k - 1] + (1.0 - alpha[i][k]) * gamma[i - 1][k - 1];

            QPointF p1, p2;

            p1 = W[i    ][k - 1] * (alpha[i][k] * gamma[i][k - 1] / gamma[i][k]);
            p2 = W[i - 1][k - 1] * ((1.0 - alpha[i][k]) * gamma[i - 1][k - 1] / gamma[i][k]);

            W[i][k] = p1 + p2;
        }
    }

    return W[j + degree][degree];
}

template< class IPoint >
vector< IPoint > NURBS< IPoint >::ComputePoints( int n ) const
{
    if ( this -> control_points.size( ) <= static_cast< unsigned int >( degree ) )
    {
        vector< IPoint > pts( n, IPoint( 0.0, 0.0 ) );
        return pts;
    }

    //int num = ( num_of_knots - 1 ) * n;
    vector< IPoint > points( n );

    for ( int i = 0; i < n; i++ )
        points[i] = Eval( static_cast< double >( i ) / static_cast< double >( n - 1 ) );

    return points;

}

template< class IPoint >
void NURBS< IPoint >::UpdateKnots( )
{
    for ( int i = 0; i <= degree; i++ )
        knots[i] = 0.0;

    for ( int i = degree + 1; i < num_of_knots + degree - 1; i++ )
        knots[i] = static_cast< double >( i - degree ) / static_cast< double >( num_of_knots - 1 );

    for ( int i = num_of_knots + degree - 1; i < num_of_knots + 2 * degree; i++ )
        knots[i] = 1.0;
}

template< class IPoint >
int NURBS< IPoint >::FindIndex( double t ) const
{
    for ( int i = degree + 1; i <= num_of_knots + degree - 1; i++ )
    {
        if ( t <= (knots[i]) )
            return i - degree - 1;
    }

    return -1;

}

#endif // NURBS_H
